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ABSTRACT 

We present a simple and efficient method to set up spherical structure models for TV- 
body simulations with a multimass technique. This technique reduces by a substantial 
factor the computer run time needed in order to resolve a given scale as compared to 
single-mass models. It therefore allows to resolve smaller scales in iV-body simulations 
for a given computer run time. Here, we present several models with an effective 
resolution of up to 1.68 x 10 9 particles within their virial radius which are stable 
over cosmolo gically re l evant time-scales. As an application, we confirm the theoretical 
prediction by lDehnenl (|2005h that in mergers of collisonless structures like dark matter 
haloes always the cusp of the steepest progenitor is preserved. We model each merger 
progenitor with an effective number of particles of approximately 10 8 particles. We 
also find that in a core-core merger the central density approximately doubles whereas 
in the cusp-cusp case the central density only increases by approximately 50%. This 
may suggest that the central region of flat structures are better protected and get less 
energy input through the merger process. 

Key words: methods: TV-body simulations - methods: numerical 



1 INTRODUCTION AND MOTIVATION 

Resolution is a key issue in TV-body simulations. In state-of- 
the-art cosmological iV-body simulations today, structures 
can be fully resolved dow n to the scale of a fract ion of a per 
cent of the virial radius l|Diemand et al. l l2007al V But there 
are many problems where even higher resolution in central 
regions of structures is needed. Also, one often would like to 
study a certain problem without the cosmological context, 
i.e. one needs a possibility to set up isolated structures with 
high resolution in order to study dynamical effects in detail 
and with precise control of the initial condition, which can 
be very difficult within the framework of a cosmological N- 
body simulation. 

For example, the question of whether the central dark 
matter density profile of haloes that form in cosmological 
iV-body simulations is cuspy or cored needs at least a res- 
olution of ss 10 -3 r v ir to be answered. Another example 
are centrally flat profiles: in order to resolve isolated struc- 
tures with flat central profiles, a lot of particles are needed 
since in flat profiles the resolution scaling with the number 
of particles is the slowest (see below for more details about 
the scaling of resolution with the number of particles) . An- 
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other problem we had in mind when developing the multi- 
mass technique presented in this paper, was the dynamics of 
super-massive black hole binaries in the centre of remnants 
of galaxy mergers. There, the relevant scales are of order of 
a few pc ~ 1CP 6 — 10~ 5 r v i r for Milky Way size dark matter 
haloes. 

We illustrate the resolution problem in more detail with 
a commonly used family of spherically sy mmetric de nsity 
profiles: the so called a/37-models family |Zhao 1996]). An 
a/37-model density profile is given by 

Po 



p(r) = 



&) 1+ & 



(1) 



where 7 determines the inner slope and /3 the outer slope 
of the density profile, whereas a controls the transition be- 
tween the inner and outer region. The normalisation is given 
by po and r s is the scale radius defined by r s = r v j r /c v ir (c v ir 
is the virial concentration) . Many models used by observers 
as well as theorists to describe structures in th e universe 
belon g to this family, e. g. the Hernquist pro file (|Hernquistj 
Il99d b the NFW profile dNavarro et all 1 19961 ) or the Moore 
profile (|Moore et al.lll998f ) are just special cases of the gen- 
eral form ([T]l. 

An obvious minimal criterion for a given length scale 
to be resolved is, that it has to be populated with enough 
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Figure 1. Plot of n, rioo, 7*relax(idyn(?*vir)/10), r re i ax (<dy n ('"vir )) as a function of iV v i r for a spherical structure with c v i r = 10, a = 1, 
/3 = 3 but different central profiles: 7 = 1 (left) and 7 = 3/2 (right). One can nicely see that in general the constraint on r TeB is set by the 
scale where relaxation becomes important since r res = max(rioo (-A^vir) i '"relax (^0>JV'vir)) f° r a given simulation time to- The asymptotic 
scaling for rjv and r re i ax is given in the plots; the discrepancy between the sampled scale rioo and the relaxation scale r^iax as a function 
of TV v ; r is bigger for larger 7. 



particles by the sampling procedure. For example, if one 
would like to set up an isolated structure given by the density 
profile {TJ, then the radius rjv that contains TV particles is 
simply given by the solution of 



M(r N ) = Nm , 

where M(r) is the enclosed mass given by 



M(r) 



4:Tvp(x)x dx , 



(2) 



(3) 



and m is the mass of the particles which can be expressed in 
virial quantities for single-mass models as m = M v i r /TV v ir- 
In the central asymptotic regime (i.e. rjv <C r s ), equation 
((2| can be solved for rjv to yield 



r N = 



3 — 7 Nm 
Air p a r] 



(4) 



For N = 1 one obtains the scale of the location of the inner- 
most particle, n mp = riQ But of course one particle is not 
enough to resolve that scale. If we say at least 100 particles 
are needed in the innermost bin so that the bin is well re- 
solved, then rioo provides a better estimate of the resolved 
scale. 

Another constraint to resolution comes from the fact 
that often the structures that one would like to simulate 

1 An alternative definition for the location of the innermost par- 
ticle can be made with the mean particle separation given by 
h(r) = The two definitions are essentially equivalent and 



only differ by the geometric factor 
unity for most useful profiles. 



477 

3-7 



which is of order 



can be very well approximated as collisionless systems, i.e. 
the local relaxation time-scale is much longer than the 
age of the system. This property must be preserved in TV- 
body simulations of such collisionless systems. Otherwise, 
purely artificial relaxation processes due to under-resolving 
the structures will have a dominant effect on the dynam- 
ics of the system dPower et al.l l 2003l ; iDiemand et al. [ |2004l ; 
iBovlan-Kolchin fc Mai 120041 ) . This is not a numerical arte- 
fact since a real astrophysical system with a small particle 
number will also be subject to relaxation processes. The 
problem is, that in most cases, we can not simulate the as- 
trophysical system under study with the number of parti- 
cles the system would have in nature, e.g. a galaxy size dark 
matter halo would have of the order of O(10 67 ) dark mat- 
ter particles if we assume that the dark matter particle has 
a mass of 100 GeV/c 2 . Therefore, this effect due to under- 
resolving the system with not enough particles will always 
be a limitation of collisionless TV-body simulations. 

In principle, there are only few known dynamically 
stable systems in the universe, e.g. the Kepl erian 2-body 
system or the figure-8 rotation of 3 bodies (|Moorel 1 19931 ; 
[Ch enciner fc Montgome ry 2000)0 Dynamical effects like re- 
laxation or evaporation will sooner or later lead to the dis- 
ruption of any system. The issue is always on what time- 
scale the system is actually stable. The time-scale of in- 
terest is normally of the order of the age of the universe 
and for most cases collisionless astrophysical systems can 

2 This is only true in the Newtonian regime. In General Relativ- 
ity, not even a 2-body orbit is dynamically stable and the orbit 
decays due to the emission of gravitat ional radiation dDvsonll979l ; 
llslamlll977l;lAdams fc Lauehlm||l997r) . 
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be regarded as perfectly stable within that time-scale. This 
is different in iV-body simulations. Here, one just tries to 
generate models that show the desired stability behaviour 
during the time-scale of interest with the least amount of 
particles needed. Artificial iV-body models will be subject 
to such disruption effects much sooner than the real astro- 
physical system one wants to study. 

We illustrate this in more details now. The local relax- 
ation time is defined by 

N(r) 



ln(jV(r)) 



where 



tdyn(r) = 2TY 



GM{r) 



(5) 



(6) 



is the dynamical or orbital time at radius r and N(r) = 
Mir) jm denotes the number of particles within r. Here, 
we have a slightly different normalisation than in the usual 
expression for the local relaxation time since we dropped 
the factor of 8 in front of the logarithmic term in the de- 
nominator. We found better agreement with results from 
iV-body simulations with this normalisation (see section 

below for more details). Had we kept the factor of 8, 
then our definition ((5j for spherically symmetric structures 
would be equivalen t to the empirical expression found by 

1 Power et al.l ([2003) . Normally, the Coulomb logarithm is 
given by ln(A) = ln(fe max /6 min ), where 6 max and 6 min are the 
maximum and minimum impact parameters of the particles 
under consideration. Since the minimum impact parameter 
is related to the softening length and the latter scales with 
the number of particles, we prefer the direct formulation of 
the Coulomb logarithm as a function of the number of par- 
ticles, ln(jV(r)). 

In a simulation run for a time to relaxation processes 
become important on a scale r re i ax given by the solution of 



^relax (^relax) — to 



(7) 



In the central asymptotic regime (i.e. r rc i ax <C r s ), equation 
([7]) can be inverted for r rc i ax to yield 



T'relax(io) = 



\ 1 g — t n I 3-7 
\ 2 3-7 to lGw7 



(8) 



where W -i denotes the fc = — 1 branch of the Lambert W 
function (|Lambertlll758l : ICorless et a.1.11 1996h and 



X 



- 7 TV 



3-7 



(3 - 7)m\ 2 < 3 -^> 



(9) 

2 3 — 7 to \Gnparl J \ ^nporl J 

In Fig. Q] we plot n, r l0 o, ■r rc i ax (t dyn (r vir )/10) and 
JYeiajc^dyn^vir)) as a function of N v i r by setting m = 
M v i r /iV v j r . The virial radius r v i r is defined so that the en- 
closed average density within r v ir is given by 

M vir 
4^ ir /3 



= pvir = A vir p cri t,o = 1.41 x 10 4 M kpc 3 (10) 



where A vir = 178 « 104 ijEke et alJ 11998ft for our 

choice of cosmology with Qm,o = 0.3, Qa.o = 0.7, p C rit,o = 
3Hi/8nG and Ho = 70 km s" 1 Mpc" 1 [ho = 0.7), which we 
use throughout this paper. The dynamical time at the virial 
radius is then 



^dyn(^"vir) — 



3tt 



GA v irPcrit,0 



= 12.2 Gyr 



(11) 



For example for a galaxy size dark matter halo with M v j r = 
W 12 /h M Q = 1.43 x 10 12 M Q one would obtain r vir = 
289 kpc and the values for systems of different virial mass 
M v i r can be obtained by the simple scaling relation 



289 kpc 



M vir 



1.43 X 10 12 Mr; 



(12) 



; max 



V ioo (iVvir ) , r rc i ax (t , N^r)) 



Although the virial radius is a rather artificial definition 
of the size of a dark matter halo and the virialised region of 
haloes in cosmological ./V-body simulations is gen erally much 
larger l|Prada et alJ feoOG; Dic mand et al. 2007b) it is a con- 
venient normalisation and cut-off scale for isolated models 
since we are anyway mainly interested in the central dynam- 
ics of the structure. We set c v i r = 10, a = 1, /3 = 3 in both 
cases and present plots for 7 = 1 (left) respectively 7 = 3/2 
(right). 

We only plot these quantities for Avir ^ 10 6 since 

the expressions Q and |SJ are only valid in the central 

_ i 

asymptotic regime. One can see that rjv oc A vir 3 ~ T and 

2 

)"rci ax ~ N vir 6 ~ J since the Lambert W function only varies 
very slowly as a function of N v ir which reflects the weak 
dependence of the logarithmic term of the local relaxation 
time. For a given particle resolution iVvir and simulation 
time to, the radius r re a that we can still resolve with correct 
collisionless physics (i.e. this scale does not suffer from too 
much artificial relaxation) is given by 

(13) 

As one can see from Fig.[T] this resolution scale r ros is in gen- 
eral set by r roi ax (to, Avir) for isolated high resolution struc- 
tures. It is worth remarking here, that for structures assem- 
bled hierarchically in a cosmological iV-body simulation, the 
amount of relaxation is signi ficantly larger and r roi ax scales 
slower as a function of A^ir l|Diemand et al ]|2004ft . Although 
a structure might be sampled with enough particles at a 
certain scale at the final time, the relaxation time at that 
scale would have been much smaller in the past since par- 
ticles were in lower mass structures during the hierarchical 
growth. 

By inspecting Fig. [T] we see that more than approxi- 
mately of order O(10 12 ) particles in the centre of a struc- 
ture with 7 = 1 are needed in order to resolve scales of 
w 10~ 5 r v ir- It is worth remarking, that with the same num- 
ber of particles iV v ir much smaller scales are populated in 
the 7 = 3/2 profile than in the 7 = 1 profile. Generally, the 
steeper the central profile, the more the particles are con- 
centrated. But unfortunately, the relaxation scale r rc iax does 
not scale equally fast so that the discrepancy as a function 
of Avir becomes bigger for larger values of 7. 

Nonetheless, such an enormous amount of particles per 
structure is hardly doable today - even with large supercom- 
puters. But since we only need this high resolution at the 
very centre of the structure, our solution to this problem is 
to use models where we only populate regions of the phase 
space that are in the centre or will reach the centre in the 
future with high resolution particles. 

In section [2] we present the simple idea behind the mul- 
timass models and present stability tests in section [3] In 
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section [4] we test Dehnen's prediction with high resolution 
mergers and we summarize our results in section [5] 



2 METHOD 

2.1 General model characteristics 

We restrict ourself to models of the form given by equation 
{TJ with 7 < 3 for the mass not to diverge in the centre. 
Similarly, for f3 ^ 3 the total mass would diverge and we 
need to introduce a cut-off radius. We chose the form 



p(r) = 



(^) 



"•decay 



r < r cut off 



r > rcutoff 



(14) 



where S and rdecav are free parameters (jSpringel fc White] 
ll999l ; lKazantzidis et ajj [2006i By requiring the logarithmic 
slope to be continuous at r cuto ff, we get 



<5 = 



7 + fl f reams. V 

cutoff ' r V '' 3 / 



I'dc 



1 



(15) 



We set the truncation scale rd CC ay = 0.3 r C utoS in order 
not to make the truncation too sharp. A too sharp trun- 
cation ca n lead t o an i nstability of the model around r C utoff 
as seen m IZemd l|2003h . For (3 > 3 we simply set r C utoff = oo 
(i.e. no cut-off) while for j3 ^ 3 one has to specify a cut-off 
scale rcutoff (e.g. r cu toff = r V ir). By further specifying r s and 
A^(r cu toff), the normalisation po is given by 

M(r cuto ff) 



po 



(16) 



where 



r fM) r ) 
dx M^JiLJ (17) 



with q = r C utoff /r s and Y is the standard gamma function 



2.2 Distribution function as probability density 
function 

For spherical systems with an isotropic velocity distribution 
one can calculate the distribution function, which in that 
case only depends on ene rgy, by the Eddington inversion 
i|Binnev fc Tremainej Il987h . Hence, by restricting to mod- 
els with an isotropic velocity distribution, we can calculate 
the distribution function for our spherical structure models 
described by equation (|14|l - at least numerically. Since the 
state of a system at a given time is completely described by 
the distribution function f(r, v), we use it as a probability 
density function in order to sample the phase space with 
particles, 



P6T>(r, v)drdv 



A/tot 



drdv 



(18) 



is the probability that a particle is in the volume drdv 
around the phase space point (r, v). By integrating out the 
velocities and using spherical symmetry, we get for the prob- 
ability density p(r) in coordinate space 



p(r)dr 



47rr 2 p(r) 
Mtot 



dr 



(19) 



and the positions can now be sampled by using the quantile 
function, which is the inverse of the cumulative probability 
distribution function M(r)/M to t for the above probability 
density function p(r). For a particle at location r, we get 
now the following probability density for the magnitude of 
the velocities 



, x, 4-Kv 2 f(ri,v) 
p{ri,v)dv = y\ — -dv 

Pin) 



(20) 



where n = \ri\ and isotropy in velocity space was used. 
In general, the distribution function can only be calcu- 
lated numerically for the large family of models described 
by the density profile (11 41) although a few analytical so- 
lutions are known (e.g. Herncmistl Il990l ; iDehnenl 1 19931 ; 
iTremaine et al.lll994T ). Hence, numerical integration and in- 
version for p(ri, v) in order to calculate the quantile function 
is difficult and one generally u s es the acceptance-rejection 
technique l|von Neumann! [19511 ; iKuiiken fc Dubinski| [l994') 
for the Monte Carlo sampling of the velocities. 

This Monte Carlo sampling procedure directly from the 
distribution function f(r, v) le ads to perfectly stable e qui- 
librium models as was shown in lKazantzidis et al.l (|2004l ) for 
single-mass models. These models do not show the flattening 
in the central part of the structure during evolution as it is 
obtained in the case of the assumption of a local Maxwellian 
velocity distribution with the velocity dispersion given by 
the Jeans equation. 



2.3 Multimass refinement 

Since in general we only need the high resolution sampling 
for structures in the central region, we use a multimass tech- 
nique to accomplish that. The general idea is that the cen- 
tral region is populated by lighter particles and the outer 
parts of the structure are sampled by heavier particles. In 
addition, a refinement depending on the orbit of the particle 
is applied in order to minimise the perturbations of heavy 
particles from the outer parts of the structure to the high 
resolution centre. 



2.3.1 Shell refinement 

By specifying an inner shell radius r s i and a number of par- 
ticles No within that radius, the mass of the particles in the 
central sphere is simply mo = M(r B i)/No and sets the ef- 
fective central resolution of the structure. One can further 
choose an outer shell radius r so , the number of shells A^hetf 
between r s i and r so , and the mass ratio between the particle 
mass in neighbouring shells A_Rm- Since one has to split par- 
ticles if one uses the orbit dependent refinement (see below 
for more details) and we do not want to have particles with 
masses smaller than mo, A_Rm has to be a natural number, 
i.e. ARm G N. The resulting structure (without orbit depen- 
dent refinement) has then a total of iV s hcn + 2 different mass 
species with rrii — mo(Ai?M)\ i = . . . iV s hcn + 1. The shells 
are equally spaced in logarithmic intervals which determines 
the number of particles within each shell. 

Of course this technique introduces some further nu- 
merical artefacts as e.g. heating of the light particles by the 
heavy particles or mass segregation of the heavy outer par- 
ticles during time evolution. The dynamical friction force 
that a particle of mass M experiences in a homogeneous sea 
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of light particles with mass m < M is Fa oc M 2 and the 
time-scale for this p article to reach the centre of the struc- 
ture is tes oc M~ _1 (jBinnev fc Tremaind Il987h . Hence, by 
choosing moderate mass ratios ARm between neighbouring 
shells, we can reduce the effect of mass segregation. By scal- 
ing the softening lengths of the heavy particles as a function 
of their mass, artificial 2-body scattering can be reduced. 
The actual role of the softening length is a cut-off scale for 
the singularities intr oduced by the Monte Carlo sampling 
jLeeuwin et alj|l993h . A natural inner scale for a given den- 
sity profile and mass resolution is set by the location of the 
innermost particle n mp = r\ oc m 3 — > (see equation Q). 
Hence, we chose to scale the softening like ri(m;) for each 
species, resulting in 



of the high resolution central region by heavy particles from 
the outer region, e.g. a heavy particle on a perfectly radial 
orbit would be split up into copies that have the same mass 
as the particles in the central high resolution region. There- 
fore, this procedure generates a protection layer around the 
high resolution region and very few heavy particles will dif- 
fuse into that region. 

2.3.3 Choice of parameters - some guidelines 

We can express the mass of a particle in shell with index i 
before the orbit dependent refinement as a function of radius 



C() 



(21) 



This can be seen as a generalised scaling for arbitrary cen- 
tral profiles 7 in isolated structures of the often used rule 
for cosmological iV-body simulations, where one scales the 
softening length by e» oc $/rrH. Of course, for 7 = 0, which 
would correspond to the homogeneous case, the two scaling 
relations are identical. 



2.3.2 Orbit dependent refinement 

In addition, we refine particles in the outer parts of the struc- 
ture depending on their orbit. From the initial position and 
velocity, we can calculate the pericentre distance r per i,fc of a 
specific particle in the smooth potential given by the density 
profile. By choosing the maximum orbital refinement radius 
r mol we determine the split factor fk for that particle by 



m 
1 

m k 1 



m 



1 °e(''pcri,fc/ r si) 

log(r mor /r si ) 
lo e( r pcri,fc/ r 3i) 

log(r j / r B i) 



^*peri,fc ^ ^si 
fmor ^ ^peri,fc 



^ _ m k 
1 



m k \ 

mo J 



< r 

< r 



peri.fe < ' mor 
pcri,fc ^ Ti ^ ^mor 



(22) 



where r; is the outer boundary of shell i that contains the 
particle with index k and is given by 



(23) 



for ^ i ^ iV s heii and rjv shcll +i = 00. If this factor fk > 1, 
then we replace that particle with f sp ut t h = int(/fc) copies 
of mass rriapiit,*; = mk/ fs P ut,k randomly placed on a sphere 
of radius dk, which is the distance of the original particle 
to the geometric centre. Here, int(:r) is the function that 
rounds the real value x to its nearest integer. We split the 
velocity of the original particle into a radial and tangential 
part. The new particles will have the same radial velocity 
component as the original particle but a new, random tan- 
gential component of the same magnitude as the original 
one. With such a splitting procedure, we keep the velocity 
configuration of the structure, i.e. the total kinetic energy is 
the same and the new particles are on the same orbit as the 
original particle. Of course, we also scale the softening of the 
new particles according to equation (|2ip . With this orbit de- 
pending refinement technique we minimise the perturbation 



m 



where ri is the outer boundary of shell % and 



log (AJ4 



(24) 



(25) 



log(r so /r si ) 

The mass in each shell Mi scales like the enclosed mass in 
the central asymptotic regime, i.e. Mi oc r 3-7 . Therefore, a 
choice of re — 3 — 7 results in an equal number of particles in 
each shell before orbit dependent refinement. Equation (JS| 

6 — ~f 

suggests, that one should scale m ~ r 2 in order to keep 
the local relaxation time constant, which would be the most 
efficient distribution of particles. But a steeper scaling with 
re > 3 — 7 also means that shells in the outer part might 
have just a few very heavy particles which does not give 
a good sampling of the structure. This behaviour becomes 
worse since the asymptotic scaling of Mi oc r 3 " 7 is only valid 
in the asymptotic central regime and becomes flatter in the 
outer parts of the density profile. Additionally, most of the 
work load is soon dominated by the high resolution centre, 
so that even fewer particles in the outer part does not result 
in a significant computer run time gain (see e.g. Fig- . 

The value of re is degenerate in the sense that many dif- 
ferent choices of iV s h c ii, r s i, r BO and A_Rm can give the same 
value for re. Here we give some guidelines on how to chose 
the different parameters. Unfortunately, there is not a set 
of universal parameters that works for all possible models 
within the a/37-family. It often also depends on the specific 
needs of the simulation. In general, we prefer not too big 
mass ratios between neighbouring shells in order to keep 
mixing and mass segregation effects to a minimum. This 
means that in most cases we set A_Rm = 2. The inner shell 
radius r s i is set so that the central region of interest is well 
sampled with high resolution particles. The outer shell ra- 
dius r so should not be chosen too large, i.e. on scales where 
one is far from the central asymptotic scaling of the enclosed 
mass. We made good experience with choosing r so ^ r v i r . 
With the final degree of freedom, the A^heii parameter, we 
control re so that re is around 3 — 7. In order to protect 
the central high resolution region from heavy particles, one 
has to chose r mor large enough. We recommend a value of 
r mol > 10 r s [. The chosen values of the stable test models 
in section [3] can also help to guide the reader choosing the 
parameters. 

Hence, with a careful choice of parameters, perturba- 
tions of structures can be minimised and effects on global 
characteristics like the radial density profile are small. The 
multimass technique is therefore an efficient method to per- 
form high resolution iV-body simulations. 
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10" 3 10" 2 10" 1 10° 10" 3 10" 2 10" 1 10° 



r / r . r / r . 

I vir I vir 

Figure 2. Multimass stability tests for mass profiles with different inner slope 7 = 0... 3/2. Each run was evolved for 10 Gyr. The top 
panels for each slope 7 show the total density profile and the sub-profiles [(light, heavy) = (dotted, dashed)] for the different mass ratio 
runs normalised to the virial density p v ; r = 1.41 X 10 4 Mq kpc - 3 . The lower panel shows the relative change of the total profile with 
radius (p — Pic)/pic normalised to the profile of the initial conditions pjc- The total density profiles are stable for moderate mass ratios 
A/?m and only show very small perturbations. For too high mass ratios, the total density profile is stronger perturbed by heating effects 
of the heavy particles and the heavy particles sink to the centre due to the higher efficiency of dynamical friction. 



Table 1. Summary of parameters for the different two-shell models. 



7 





1/2 


1 


3/2 


eo [r V ir] 


3.11 x 10~ 3 


1.56 x 10~ 3 


1.04 x 10~ 3 


5.18 x 10~ 4 


N ■ 

1 v vir 


7.21 x 10 6 


4.88 x 10 6 


3.25 x 10 6 


2.11 x 10 6 


»"rclax(10 Gyr, 7V vir ) [r v ir] 


2.44 x KT 3 


2.25 x 10~ 3 


2.02 x 10~ 3 


1.78 x 10~ 3 
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3 TESTS 

3.1 Two-shell models without orbit dependent 
refinement 

In order to illustrate the basic idea and the limitations of the 
multimass technique, we present a series of simple two-shell 
models without an orbit dependent refinement, i.e. r mor — 0. 
In this series we chose four different profiles from the aj3^- 
family described by equation (|14|l . The following density 
profile parameters were the same for all models: outer pro- 
file (5 = 3, transition coefficient a = 1 and concentration 
c v i r = 20. We varied the inner profile form 7 = ... 3/2. For 
the refinement we chose r s i = r so = r s , iV s heil = (i.e. only 
two shells) and iVo = 3 x 10 5 for all models. For the different 
runs, we varied the mass ratio ARm with the following val- 
ues lx 10°, 3x 10°, lx 10\ 3x 10\ Ixl0 2 ,3xl0 2 , lxlO 3 . The 
choice of softening for the high resolution particles eo and the 
total number of particles jV v ir in the case of equal mass par- 
ticles in the inner and outer shell, which corresponds to the 
effective resolution of the multimass models, as well as the 
estimated relaxation scale after 10 Gyr, r rG i ax (10 Gyr, N v i T ), 
are given in Table [T] The softening lengths of the heavy 
particles were scaled as described by equation (|21|) . 

Each of these 28 models was evolved in isolation for 
10 Gyr in order to test the stability of the structures. For 
the time evolution of all models we present in this paper, 
we use pkdgrav, a state-of-the -art tree code written by 
by Joachim Stadel ( Stadcl 2 00 ll ) . We used the dyna mical 
time-stepping scheme developed bv lZemp et al.l l|2007h . For 
multimass models it is important to use a time-step for the 
particles that is based on the true dynamical time of a par- 
ticle and independent of the softening because heavy parti- 
cles with a large softening and light particles with a small 
softening can now mix. If for example two such particles 
were at the same geometric distance from the centre they 
should both take the same time-step. This is guaranteed by 
a criterion that is based on the dynamical time. Time-step 
criteria that are based on or scale with characteristics of the 
particle can't fulfil that and can lead to physically wrong 
time-steps. In addition, for high resolution models, the dy- 
namical time-stepping is more efficient and more accurate in 
the high resolution centr e. For a detailed d iscussion about 
time-stepping criteria see IZemp et all (|2007T l. 

In Fig. [21 we present the density profiles of these runs 
after 10 Gyr. For moderate mass ratios ARm up to 10-30 
(or even ARm ~ 100 for steep central profiles) the effects 
on the total density profile are small and the profile remains 
stable down to the level of a few eo. Such deviations are any- 
way expected since the forces in pkdgrav are softened if two 
particles have distances of order of their softening length. By 
comparing the equal mass cases for the two steepest central 
profiles, we see that the flattening effect due to relaxation 
sets in at a radius that is a little bit smaller than the esti- 
mated value r ro i ax (10 Gyr, iV v ir)- Hence, our estimate ((5J) is 
a good estimate. If we would keep the factor 8 in front of 
the logarithmic term for the relaxation scale, r ro i ax would be 
approximately twice as large, confirming that it is rather a 
conservative estimate. This was our motivation to drop the 
factor of 8. 

Figure [5] also illustrates that the different species form 
stable sub-profiles. The shell radius was chosen in a zone 
where the local density profile is steep so that the transition 
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Figure 3. The theoretical and measured speed up factors Sth 
respectivley s m as a function of mass ratio Ai?M for the different 
central profiles 7 for the two-shell models without orbit dependent 
refinement. A substantial fraction of time is gained in all cases. 
Most of the time gain is already obtained for small mass ratios. 



region is small and in the inner or outer regions the light 
or heavy particles dominate respectively. Only for very high 
mass ratios, the total density profile is strongly perturbed by 
heating effects of the heavy particles and the heavy particles 
sink to the centre due to the higher efficiency of dynamical 
friction. But for moderate mass ratios, these effects are small 
and expected only for much longer time-scales. 

The main advantage of these multimass models is the 
speed up gain. We can estimate a theoretical speed up factor 
Sth by the ratio of the number of force evaluations between 
a single-mass and a multimass model of the same density 
structure 



x> 47rp(r)H , 
imp m AT(r) U ' 



Z^i=0 Jr i _ 1 

E 



47rp(r)7' 2 
iHjAT(r) 



dr 



j=l m AT(rj) 



AT(rj) 



where r» is the outer boundary of shell 



and we set 



and 



(26) 
(27) 

(see also equation 
1 = 00. AT(r) is 



the time-step a particle takes at radius r. Here we use a 
dynamical time-stepping scheme where 



AT(r) = 7? D 



r- 5 _ no_ , . 
GM{r) ~ 2tt dynlrj 



(28) 



with ?7d = 0.03. For more detai ls about the dynam ical time- 
stepping scheme please consult TZemp et al.l (2007). The sec- 
ond expression for sth is valid for a discrete sampling of the 
structure and can be obtained from the first expression by 
plugging in the discrete density 
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KO = f>^ (29) 

3=1 

where we used the Dirac delta distribution in spherical co- 
ordinates. 

In Fig. |31 we plot the theoretical estimate s t h an d the 
measured speed up factor s m = Tq/T, where To is the time 
needed by the equal mass run for the 10 Gyr simulation time 
stability test run described above and T is the computer run 
time needed by the multimass model. Figure [3] illustrates 
that we gain a substantial fraction of computer run time in 
all cases. Most of the gain is already obtained for small mass 
ratios Rm- For example, the run with inner slope 7 = and 
mass ratio ARm ~ 10 is approximately four times faster 
than the same run without multimass refinement and does 
not show any perturbation effect of the multimass technique 
on the density profile. 

The steeper the central profile, the less is the computer 
run time gain. This is due to the fact that most of the work 
in the iV-body simulation for steep profiles is concentrated 
in the centre. In the centre, the particles are on very small 
time-steps compared to the less dense, outer regions of a 
structure in an iV-body simulation and as a consequence a 
lot of expensive force calculations are needed. 

For flat density profiles the theoretical speed up esti- 
mate Sth agrees quite well with the measured value s m . Only 
for steeper central profiles the measured speed up factor s m 
is larger than the theoretical prediction. This can be ex- 
plained by the fact that in the steep profiles the particles 
in the very centre are expelled due to numerical effects (see 
for example the central density profiles in figure [2}. As men- 
tioned above, in steep central profiles the work is dominated 
by the central particles. Hence, if one looses some of them by 
numerical effects the work load decreases and the speed up 
factor becomes larger than the theoretical estimate which 
assumes that the density profile stays perfectly stable over 
time. In flatter central profiles this numerical effect is much 
weaker and the theoretical estimate agrees quite well with 
the measured speed up factor. In that sense, the theoretical 
estimate s t h is a lower limit for the speed up gain. 

3.2 Multimass models with orbit dependent 
refinement 

In the previous section, we presented models where the the 
heavier particles were allowed to penetrate the centre. Now, 
we focus on more general models with several particle species 
that have in addition the orbit dependent refinement de- 
scribed earlier. We performed a series of runs for a structure 
model with a = 1, /3 = 3, 7 = 1, and c v i r = 10. All models 
had r si = 3.46 x 10 -3 r vir , iV = 10 4 , e = 2.59 x 10" 4 r vir 
in common. The other refinement parameters are given in 
Tab. [2] This results in a structure with an effective reso- 
lution within the virial radius of N^ T — 2.61 x 10 7 parti- 
cles. For comparison, we also set up a single-mass model 
with the same resolution but without refinement as a ref- 
erence (model R). Model R was sampled with a total of 
-/Vaamp = 3.54 x 10 7 particles including the cut-off region. 
We also kept A_Rm = 2 since we now have the parameter 
iVsheii to control the overall mass range of the species and 
with that choice the mass contrast between neighbouring 
shells is kept to a minimum. 



Each of these models was evolved for 10 Gyr in order 
to test for stability. The relaxation scale for such a halo is 
r rc iax(10 Gyr, N^f T ) = 1.07 x 10 -3 r vir . We present stability 
plots for the more progressive models of series B and C in 
Fig. [4] We see that in general the profiles stay stable over 10 
Gyr down to the relaxation scale and the perturbations are 
on a few per cent level - even for the reference single-mass 
model. Only for model C3 we pushed the refinement param- 
eters too far (i.e. too small values for r ao and r m0 r) and the 
structure is relatively strongly perturbed. This was done de- 
liberately in order to demonstrate that this technique has its 
limits and needs some careful parameter choice by the sim- 
ulator that generates the initial conditions. We recommend 
in general to choose the values for r so and r mor not too small 
with respect to r B i. This guarantees a larger protection zone 
around the high resolution centre and minimises the pertur- 
bations by the heavier particles. 

The main advantage is again the computer run time 
gain. As we can see in Fig. [4] many models which are stable 
down to a similar scale as the single-mass reference model R, 
are approximately an order of magnitude faster. In Tab. [5] 
we give for each model the values for Sth which we estimated 
from the discrete sampling of each structure (see equation 
and s m = Tq/T where To is the computer run time 
for the single-mass reference model R and T is the computer 
run time of the specific multimass model. We also see that a 
more progressive choice of k as for example between model 
Bl and B2 can lead to a 20% speed up. Of course, we ran 
all models under the same conditions, i.e the same number 
of CPUs on the same super-computer. 

We now illustrate the mixing of the different particle 
species in more detail on the basis of model Bl and B3. For 
Bl we chose iV s h e ii = 5 and r mor = 3.46 x 10 -2 r vlr whereas 
for B3 we chose some more progressive values of N a ^ e \\ = 10 
and r mor = 1.73 x 10~ 2 r vlI . The model B3 is then sampled 
by less than a third of the particles needed for Bl and the 
computer run time for model B3 is approximately half the 
value of Bl. If we say we trust the dynamics down to the 
scale where the relative perturbations (p — pic) /pic are not 
larger than 10%, then we lose approximately a factor of two 
in resolution for model B3 compared with Bl (see also Fig. 

In Fig. O we plot histograms of the particle distribution 
for both models, Bl (top) and B3 (bottom), for the initial 
conditions (left) and after 10 Gyr (right). The JVgheii shells 
with index 1 . . . A r s i lc n are defined by the values of r s i and 
r so . Between these two values, the shells are equally spaced 
logarithmically. The shell with index contains everything 
smaller than r s i and the shell with index iV s heii + 1 includes 
everything larger than r ao . The ranges for the mass species 
are defined in a similar way: the mass species with index i 
contains the range m;_i < m ^ mi where m; = mo(A7?M) 1 
and m_i = 0. We see that even after 10 Gyr only a few 

3 It is in principle also possible to estimate s t h in the case with 
orbit dependent refinement from a similar continuous estimate 
like equation J26D for the case without orbit dependent refine- 
ment. This estimate corrects in addition for the mass fraction at 
a given radius r that is splitted according to the splitting proce- 
dure described in section [27372] Unfortunately, such an estimate is 
rather complicated to evaluate analytically and an estimate from 
the discrete sampling of the structure is much more practical. 
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Figure 4. Stability plots for the more progressive models of series B (left) and C (right) after 10 Gyr. One can see that in general the 
profiles stay stable over 10 Gyr down to the relaxation scale and the perturbations are on a few per cent level - even for the reference 
single-mass model. Only for model C3 we pushed the refinement parameters too far and the structure is relatively strongly perturbed. 



Table 2. Summary of parameters for the different models with orbit dependent refinement. Also given are the theoretical and measured 
speed up factors for each model. 



Model 


Al 


A2 


A3 


Bl 


B2 


B3 


CI 


C2 


C3 


fso [?™vir] 


1 


1 


1 


3.46 x lO" 1 


3.46 x lO" 1 


3.46 x 10" 1 


1.04 x 10" 1 


1.04 x lO" 1 


1.04 x 10" 1 


fmor [fvir] 


3.46 x 10" 2 


3.46 x 10~ 2 


1.73 x 10~ 2 


3.46 x 10" 2 


3.46 x 10~ 2 


1.73 x 10" 2 


3.46 x 10" 2 


3.46 x 10~ 2 


1.73 x 10~ 2 




5 


10 


10 


5 


10 


10 


5 


10 


10 


N samp 


2.53 x 10 6 


1.13 x 10 6 


6.67 x 10 s 


1.95 x 10 6 


1.06 x 10 6 


6.02 x 10 5 


1.63 x 10 6 


1.05 x 10 6 


5.51 x 10 5 


K 


0.612 


1.22 


1.22 


0.753 


1.51 


1.51 


1.02 


2.04 


2.04 


8th 


4.67 


6.92 


9.76 


5.40 


7.32 


10.7 


6.30 


7.69 


11.7 


Sxn 


5.75 


8.82 


12.3 


6.61 


8.43 


13.3 


7.51 


8.72 


12.9 



particles with masses heavier than mo are in the innermost 
shell. For example in model Bl, only four particles of mass 
mi are in the innermost shell after 10 Gyr. A few more parti- 
cles from neighbouring shells populate the innermost shell in 
model B3, which has a less conservative choice of refinement 
parameters. For model Bl we chose r mor = 10 kpc which lies 
in the middle of shell 3 and for B3 we chose r mor = 5 kpc 
which lies in the middle of shell 4. We see that none of the 
very heavy particles are closer than r mor in both cases. This 
tells us that the orbit dependent refinement procedure works 
very well and we do not get scattering of heavy particles on 
orbits that go through the centre. 



3.3 A structure with 1.68 x 10 9 particles 

As a last test we present a structure with an effective res- 
olution of Ny® — 1.68 x 10 9 particles. We chose a = 1, 



(5 = 3, 7 = 1, Cvii 
r si = 5 x 10" 3 r 
k = 1.15 and r n 



= 20, N = 10 4 , e = 1.73 x 10" 



= 3.46 x 10" 
6.91 x 10~ 3 



ATshell 

. The theoretical 



= 12, 



speed up factore for this structure is Sth = 28.6 which 
we again estimated from the sampled structure via equa- 
tion (|27[) . We evolved this structure for 5 Gyr and the 
density plot can be seen in Fig. [6] The relaxation scale is 
r re iajt(5 Gyr,A^S) = 1.31 x 10 -4 r vl[ and we see that this 
structure is stable down to that scale. 



4 PRESERVATION OF CUSP SLOPE 



iDehnenl (2005) showed analytically that in a merger of self- 
gravitating cusps with different central slopes the merger 
remnant has always the slope of the steepest progenitor. In 
other words the steepest cusp slope is preserved. 

Collisionless mergers of dark matter halos 
were already studied in e arlier work (e.g. iBarnesI 



19991: iBovlan-Kolchin fc Mai 120041: iMoore etall 120041 : 



Aceves fc Velazquez! 120061 ; iKazantzidis et alj 120061 ) but 
none of these studies had the resolution of the simulations 
presented here. We do not find any discrepancies between 
this work and the earlier studies lending support to the 
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Figure 5. Histograms of the particle distribution for the models Bl 
Gyr (right). We see that even after 10 Gyr only a few particles with 
the very heavy particles are closer than r mor in both cases. 

fact that the multimass technique works under the extreme 
dynamics of a dark matter halo merger event. 

We initially set-up two dark matter halo models with 
different central slope in isolation as in section [3] Both 
profiles were from the a/37-model family and had the fol- 
lowing general specifications: a = 1, P = 3, c v i r = 10, 
M vir = 1.43 x 10 12 M Q and r vir = 289 kpc. One struc- 
ture had a central slope of 7 = 0, the other had 7=1. For 
both models we used A_Rm = 2 and the other parameters 
are given in Tab. [3] The resolution is clearly limited by the 
centrally flat (7 = 0) model since here many more parti- 
cles are needed to resolve a given scale compared to steeper 
profiles. To resolve the same scale with the steeper profile 
then fewer particles would be needed, but since we would 
not like to have too high mass ratios between the two high 
resolution species in the centre we increased the number of 
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(top) and B3 (bottom) for the initial conditions (left) and after 10 
masses heavier than mo are in the innermost shell and that none of 



particles in the 7 = 1 structure. The theoretical speed up 
factor is s t h = 12.0 for the model with 7 = and s t h = 7.39 
for the model with 7 = 1 which was again estimated from 
the discrete sampling of the structures. 

We evolved both halos in isolation for 10 Gyr in or- 
der to test the stability again. In total three mergers were 
performed: a cusp-cusp, a core-core and a cusp-core merger 
where cusp means 7 = 1 and core is equivalent to the 7 = 0. 
The following merger set-up was used for all three cases. 
We placed the two individual haloes 600 kpc ~ 2 r v i r apart 
and the two haloes had an initial relative radial velocity of 
«rad = 150 kpc Gyr -1 and a relative tangential velocity of 
vtan = 50 kpc Gyr -1 . Assuming the two haloes were point 
masses, this set-up corresponds to an eccentricity e ~ 0.95 of 
the orbit consistent with values in co smological iV-body sim- 
ulations (jKhochfar fc Burkertl 120061 ). With this set-up the 
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Figure 7. Equal slope mergers, left column: 7 = and right column: 7 = 1. Top row: The profiles are stable over fO Gyr and the central 
slope of the structure is conserved in both cases. Interestingly, we find that in the 7 = case the central density approximately doubles 
whereas in the 7 = 1 case the central density in the merger remnant only increases by approximately 50%. This can also be seen in the 
average sub-profile of the two individual structures. Bottom row: histograms of the particle distribution after 10 Gyr. We can again see 
that only a few heavy particles are in the innermost shell. 



merger time needed by the two haloes to merge completely 
was approximately 7 Gyr. We let all three runs evolve to 10 
Gyr so that the merger remnant has time to relax. 

Fig. [7] the top row shows the stability of both models 
and the equal profile mergers after 10 Gyr. Down to a few 
softening lengths of the lightest particles no significant de- 
viations can be seen and the merger remnant has the same 
central profile as the two progenitors. We normalise the plots 
by the values r v i r and p v ir of the progenitor haloes. Inter- 
estingly, we find that in the 7 = case the central density 
approximately doubles whereas in the 7 = 1 case the central 



density in the merger remnant only increases by approxi- 
mately 50%. This can also be seen in the average sub-profile 
of the two individual structures which stays constant in the 
7 = case whereas it decreases by approximately 25% in 
the cusp case with 7=1. The case with 7 = 1 is similar 
to a major merger in a cosmological simulation where one 
finds the same effect (fFa ltenbach er et al.ll2006h . This may 
suggest that the central region of flat structures are better 
protected and get less energy input from the merger. If this 
is a general feature of cored profiles and to what degree the 
merger configuration plays a role is unclear and needs to be 
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Figure 6. Density profile of an equilibrium structure with APS = 
1.68 X 10 9 particles after 5 Gyr. 



Table 3. Summary of parameters for the two models used for the 
mergers. 
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1 




N 


1 x 10 4 




4 x 10 4 




Kir] 


3.46 x 10" 


4 


3.46 x 10" 


4 


N cS 
vir 


1.14 x 10 


i 


1.04 x 10 


t 


N samp 


5.67 x 10 




4.29 x 10 




^*si L^Vir] 


6.91 x 10" 


3 


3.46 x 10" 


3 


^so [^*vir] 


3.46 x 10" 


1 


3.46 x 10" 


1 


^mor [^vir] 


6.91 x 10" 


2 


3.46 x 10" 


2 


A/shell 


10 




10 




K 


1.77 




1.51 




r ro la X (10 Gyr) [r vir ] 


1.19 x 10" 


3 


6.26 x 10" 


4 



investigated further. The bottom row shows histograms of 
the particle distribution after 10 Gyr. The particle species 
were defined in the same way as for Fig. [5] Even under the 
violent dynamcis of a merger, the contamination of the in- 
nermost shell with heavier particles is relatively small and 
hence does not affect the central slopes of the density pro- 
files. 

In Fig. [8] we present the profile of the cusp-core merger 
after 10 Gyr. The central slope of the steepest progenitor 
is perfectly preserved and the cored progenitor only con- 
tributes significantly in the outer region of the density pro- 
file. If we look at the sub-profile of the particles that be- 
longed initially to the cored halo, we now see that in this 
case the cored structure gets much more energy input from 
the merger and its central density is lowered approximately 
by a factor of three through the merger process. 

We can therefore confirm the earlier findings that core- 
core mergers lead to a cored merger remnant while cups-cusp 
mergers lead to a cuspy merger remnant with high resolu- 
tion multimass iV-body simulations. In cusp-core mergers 



Figure 8. Cusp-core merger after 10 Gyr. The central cusp is 
preserved. As a reference we also pot the initial profiles of two 
individual structures. 



the merger remnant has a final profile corresponding to the 
steepest progenitor which is in exc e llent agreement with the 
theoretical predictions bv lDehnenl (|2005h . 



5 CONCLUSIONS 

The multimass technique for modelling haloes is a simple 
method to perform high resolution Af-body simulations. An 
early version of this technique without orbit dependent re- 
finement has already been used successfully in several appli- 
cations whic h also include cosmological structure for mation 
simulations ijDiemand et al.ll2005l : iGoerdt et alJl2006h . With 
a careful choice of parameters it is possible to gain over an 
order of magnitude in computer run time for a given reso- 
lution scale. 

As an application of this technique we confirm the ear- 
lier findings that core-core mergers lead to a cored merger 
remnant while cups-cusp mergers lead to a cuspy merger 
remnant with high resolution multimass A^-body simula- 
tions. In cusp-core merges, the merger remnant has a fi- 
nal profile corresponding to the steepest progenitor which is 
in excellent agreement with the theoretical predictions. We 
find that in the core-core case the central density approx- 
imately doubles whereas in the cusp-cusp case the central 
density in the merger remnant only increases by approxi- 
mately 50%. This may suggest that the central region of flat 
structures are better protected and get less energy input 
from the merger. 

A software tool called halogen (HALO GENerator) for 
generating multimass initial conditions is available from the 
author upon request. 
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